Los datos con los que se trabajará en este TP provienen de la 3° Encuesta Mundial de Salud Escolar (EMSE) provistos por el Ministerio de Salud link de la República Argentina. Esta encuesta trata sobre temas de salud y hábitos de las personas en la escuela secundaria que pueden impactar en su salud.
rm( list=ls() ) #remove all objects
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 541512 29 1212557 64.8 NA 669240 35.8
Vcells 1037625 8 8388608 64.0 65536 4787859 36.6
options(scipen=999)
# Carga de librerías
library(corrr)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(tidymodels)
library(rsample)
library(gridExtra)
library(knitr)
library(kableExtra)
library(GGally)
Comenzamos leyendo los datos y viendo su estructura, para esto utilizamos la función glimpse
#Importamos los datasets con los que trabajaremos
encuesta_salud_train = read.csv("/Users/dfontenla/Maestria/2022C2/EEA/practica/repo/EEA-2022/TP1/encuesta_salud_train.csv", encoding = "UTF-8") %>% mutate(id = 1:nrow(.))
encuesta_salud_test = read.csv("/Users/dfontenla/Maestria/2022C2/EEA/practica/repo/EEA-2022/TP1/encuesta_salud_test.csv", encoding = "UTF-8") %>% mutate(id = 1:nrow(.))
#encuesta_salud_train$nivel_educativo
#Observamos la estructura y variables
encuesta_salud_train %>% glimpse()
Rows: 7,024
Columns: 17
$ record <int> 502, 26488, 31473, 14154, 36578, 53730, …
$ edad <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero <chr> "Femenino", "Masculino", "Masculino", "M…
$ nivel_educativo <chr> "2do año/11vo grado nivel Polimodal o 4t…
$ altura <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ frecuencia_hambre_mensual <chr> "Rara vez", "Rara vez", "Nunca", "Nunca"…
$ dias_consumo_comida_rapida <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1…
$ edad_consumo_alcohol <chr> "14 o 15 años", "7 años o menos", "Nunca…
$ consumo_diario_alcohol <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…
$ consumo_semanal_frutas <chr> "No comí frutas durante los últimos 7 dí…
$ consumo_semanal_verdura <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_gaseosas <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_snacks <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_comida_grasa <chr> "No comí comida alta en grasa en los últ…
$ id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1…
Existen 7024 registros, compuesto de 17 datos diferentes con los que trabajaremos para entrenar diferentes modelos y poder realizar la prediccion de la variable peso. Nos quedamos con todo el dataset, como se menciona en el enunciado, corresponden a un recorte (muestra) del dataset original, luego del tratamiento de valores atípicos e ingeniería de atributos.
tabla_exploratorios = encuesta_salud_train %>%
gather(.,
key = "variables",
value = "valores") %>% # agrupamos por las variables del set
group_by(variables) %>%
summarise(valores_unicos = n_distinct(valores),
porcentaje_faltantes = sum(is.na(valores))/nrow(encuesta_salud_train)*100) %>%
arrange(desc(porcentaje_faltantes), valores_unicos) # ordenamos por porcentaje de faltantes y valores unicos
tabla_exploratorios
Observamos un dataset que ya esta curado de datos faltantes y su composición esta dada por una mayor cantidad de datos categóricos que continuos
Realizamos un primer análisis descriptivo de nuestro dataset de columnas relevantes. Para ello utilizamos la función ggpairs sobre las variables numéricas con una apertura por la variable de sexo.
#Observamos la correlacion de variables
#encuesta_salud_train_numeric = encuesta_salud_train %>% dplyr::select(where(is.numeric))
# Seleccionamos las variables de interés
encuesta_salud_train_gender_numeric = encuesta_salud_train %>%
dplyr::select(edad, genero, altura, peso,dias_consumo_comida_rapida, consumo_diario_alcohol,
dias_actividad_fisica_semanal)
encuesta_salud_train_gender_numeric %>% glimpse()
Rows: 7,024
Columns: 7
$ edad <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero <chr> "Femenino", "Masculino", "Masculino", "M…
$ altura <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ dias_consumo_comida_rapida <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1…
$ consumo_diario_alcohol <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…
encuesta_salud_train_gender_numeric %>% ggpairs(aes(color=genero), upper = list(continuous = wrap("cor", size = 3, hjust=0.8, alignPercent=0.15)), legend = 25) +
theme_bw() +
theme(axis.text.x = element_text(angle=45, vjust=0.5), legend.position = "bottom")
Realizamos el análisis solo numérico de la correlación entre variables
encuesta_salud_train %>%
correlate() %>% # convierte la matriz de corr en dataframe
shave() %>% # solo muestra información debajo de la diagonal principal
fashion() # acomoda los datos en forma tidy (por ej. redondeo de decimales)
Lo graficamos para tener las relaciones de una forma mas descriptiva
encuesta_salud_train %>%
correlate() %>%
network_plot(min_cor = 0.1)
Respecto a nuestra variable a predecir, el peso, observamos:
Para las categorías de la variable frecuencia de hambre mensual, analice gráficamente la distribución en términos de frecuencia relativa de:
Agrupamos los valores para obtener la frecuencia tanto de consumo semanal de verduras como de consumo semanal de comida grasa
#encuesta_salud_train_numeric <- data.frame(encuesta_salud_train)
consumo_verduras_freq = encuesta_salud_train %>%
group_by(frecuencia_hambre_mensual, consumo_semanal_verdura) %>%
summarise(n = n()) %>%
mutate(freq = n / sum(n))
consumo_verduras_freq %>% glimpse
Rows: 45
Columns: 4
Groups: frecuencia_hambre_mensual [6]
$ frecuencia_hambre_mensual <chr> "Algunas veces", "Algunas veces", "Algunas v…
$ consumo_semanal_verdura <chr> "1 a 3 veces durante los últimos 7 días", "1…
$ n <int> 175, 92, 79, 21, 119, 27, 3, 73, 23, 14, 5, …
$ freq <dbl> 0.297113752, 0.156196944, 0.134125637, 0.035…
consumo_semanal_comida_grasa_freq = encuesta_salud_train %>%
group_by(frecuencia_hambre_mensual, consumo_semanal_comida_grasa) %>%
summarise(n = n()) %>%
mutate(freq = n / sum(n))
consumo_semanal_comida_grasa_freq %>% glimpse
Rows: 43
Columns: 4
Groups: frecuencia_hambre_mensual [6]
$ frecuencia_hambre_mensual <chr> "Algunas veces", "Algunas veces", "Alguna…
$ consumo_semanal_comida_grasa <chr> "1 a 3 veces durante los últimos 7 días",…
$ n <int> 252, 62, 33, 18, 81, 24, 4, 115, 27, 6, 9…
$ freq <dbl> 0.427843803, 0.105263158, 0.056027165, 0.…
ggplot(data = consumo_verduras_freq, aes(x = frecuencia_hambre_mensual, y = freq, fill = consumo_semanal_verdura)) +
geom_bar(stat = "identity", position = position_dodge()) +
xlab("Frecuencia de hambre mensual") +
ylab("Consumo semanal de verduras") +
theme(axis.text.x = element_text(angle = 90)) +
coord_flip() +
theme_bw() +
labs(title = "Frecuencia de hambre mensual") +
scale_fill_brewer(palette = "Paired") +
theme(axis.text.x = element_text(size = 10)) +
theme(axis.title.x = element_text(size = 15)) +
theme(axis.title.y = element_text(size = 15)) +
theme(plot.title = element_text(size = 20))
Se puede observar como para las categorias Siempre y Casi siempre de hambre mensual, aumenta de forma significativa la cantidad de valores en consumo semanal verduras, que hacen referencia a no ingerir estos alimentos con recurrencia, y si aquellas que referencian a valores semanales (1 a 3 veces durante los ultimos 7 dias y 4 a 6 veces durante los ultimos siete dias) en lugar de valores diarios
ggplot(data = consumo_semanal_comida_grasa_freq, aes(x = frecuencia_hambre_mensual, y = freq, fill = consumo_semanal_comida_grasa)) +
geom_bar(stat = "identity", position = position_dodge()) +
xlab("Frecuencia de hambre mensual") +
ylab("Consumo de comida grasa") +
theme(axis.text.x = element_text(angle = 90)) +
coord_flip() +
theme_bw() +
labs(title = "Frecuencia de hambre mensual") +
scale_fill_brewer(palette = "Paired") +
theme(axis.text.x = element_text(size = 10)) +
theme(axis.title.x = element_text(size = 15)) +
theme(axis.title.y = element_text(size = 15)) +
theme(plot.title = element_text(size = 20))
Por otro lado, en cuanto al análisis de la relación entre el consumo de comida grasa y la frecuencia de hambre, observamos que el consumo habitual de comidad grasa, aumenta la frecuencia de hambre mensual. Buscando referencias a esta conclusión, existen estudios relacionados a la obecidad donde focalizan sobre esta relación donde los alimentos con altos contenidos grasos, dan más hambre
Se plantea que una primera alternativa para modelar el peso es: E(peso) = β0 +β1altura+β2edad+β3genero+β4diasActividadF isicaSemanal+β5consumoDiarioAlcohol
Realizamos una selección de las variables de nuestro interés para poder realizar el modelado, y observamos su estructura
#encuesta_salud_train %>% glimpse
encuesta_salud_modelo_base = encuesta_salud_train %>%
# Seleccionamos las variables de interés
dplyr::select(edad, genero, altura, peso,dias_actividad_fisica_semanal,consumo_diario_alcohol)
encuesta_salud_modelo_base %>% glimpse
Rows: 7,024
Columns: 6
$ edad <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero <chr> "Femenino", "Masculino", "Masculino", "M…
$ altura <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…
$ consumo_diario_alcohol <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …
Vemos cómo es la correlación entre las variables numéricas.
encuesta_salud_modelo_base %>% ggpairs(aes(color=genero), upper = list(continuous = wrap("cor", size = 3, hjust=0.8, alignPercent=0.15)), legend = 25) +
theme_bw() +
theme(axis.text.x = element_text(angle=45, vjust=0.5), legend.position = "bottom")
Armamos un modelo para predecir el peso en función de la altura, edad, género, días de actividad física semanal y consumo diario de alcohol.
E(peso) = β0 +β1altura+β2edad+β3genero+β4diasActividadFisicaSemanal+β5consumoDiarioAlcohol
# ajustamos modelo lineal múltiple
modelo_ln_base <- lm(peso ~ edad + genero + altura + dias_actividad_fisica_semanal + consumo_diario_alcohol, data = encuesta_salud_modelo_base)
# Resumen del modelo
tidy_ln_base <- tidy(modelo_ln_base, conf.int = TRUE)
tidy_ln_base
El valor de β0^ (ordenada al origen) es -68.92 kg, lo que corresponde al peso esperado de una persona femenina sin edad, sin altura, sin dias de actividad fisica y sin consumo diario en alcohol. Lo cual, este caso carecería de sentido ya que las personas femeninas deberían tener como valores de altura y edad numeros positivos > 0.
β0^ + βgeneroMasculino es la media del peso para las personas de género masculino, dada la edad, altura, dias de consumo de alcohol, dias de actividad fisica. Por lo tanto, βgeneroMasculino es la diferencia en los niveles medios de pesos de las personas masculinas respecto de las femeninas (categoría basal). Es decir, βgeneroMasculino (1.262643558) indica cuánto más alta es la función de respuesta (peso) para las personas masculinas respecto de las femeninas (categoría basal), dada la edad, altura, dias de consumo de alcohol y dias de actividad física.
El coeficiente estimado de edad es de 1.4kg, lo que indica que si mantenemos el número de altura, género, dias de actividad fisica y consumo diario de alchol, cada incremento adicional de edad corresponde a un aumento de 1.4kg, en promedio en el peso de la persona femenina. O lo que es igual, dadas dos personas con la misma altura, dias de actividad fisica, genero y consumo de alcohol, pero teniendo una un año más de edad que la otra, el peso esperado para la de mayor edad será 1.4kg más que la de menor pesaje.
El coeficiente estimado de altura es de 0.65kg, lo que indica que si mantenemos el número de género, edad, dias de actividad física y consumo diario de alchol, cada incremento adicional de altura corresponde a un aumento de 0.65kg, en promedio en el peso de la persona femenina. O lo que es igual, dadas dos personas con la misma edad, dias de actividad fisica, genero y consumo de alcohol, pero teniendo una un centímetro más de altura que la otra, el peso esperado para la de mayor altura será 0.65kg más que la de menor pesaje.
El coeficiente estimado de dias de actividad fisica es de -0.0874kg, lo que indica que si mantenemos el género, edad, altura y consumo diario de alchol, cada incremento de un día de actividad fisica adicional corresponde a un decremento de 0.0874kg, en promedio en el peso de la persona femenina. O lo que es igual, dadas dos personas con la misma edad, altura, género y consumo de alcohol, pero teniendo una un día de actividad fisica más que la otra, el peso esperado para la de mayor dias de actividad física será 0.0874kg menor que la de mayor pesaje.
El coeficiente estimado de dias de consumo de alcohol es de 0.00727kg, lo que indica que si mantenemos el número de género, edad, altura y días de actividad física, cada incremento de un día de consumo de alcohol adicional corresponde a un aumento de 0.00727kg, en promedio en el peso de la persona femenina. O lo que es igual, dadas dos personas con la misma edad, altura, género y días de actividad física, pero teniendo una un día de consumo de alcohol más que la otra, el peso esperado para la de mayor días de consumo de alcohol será 0.00727kg más que la de menor pesaje.
Para evaluar la significatividad individual de cada una de las variables se analiza el test t que busca probar si el coeficiente de regresión correspondiente a dicha variable es distinto de 0 (figura en la tabla resumen de resultados de la regresión).
Es decir, buscamos probar:
H0:βk=0 H1:βk≠0.
options("scipen"=1)
tidy_ln_base %>% select(term, statistic, p.value, conf.low, conf.high)
significancia = tidy_ln_base %>% dplyr::select(term, statistic, p.value, conf.low, conf.high)
significancia
Evaluando la significancia de cada variable observamos lo siguiente
Para las variables edad, altura, generoMasculino, generoFemenino (basal):
Para las variables días de actividad física semanal y consumo diario de alcohol:
tidy(anova(modelo_ln_base))
La tabla de ANOVA muestra que, según el resultado del test F, la variable género en su conjunto, también resulta estadísticamente significativa para explicar al precio (p-valor < 0.05).
Se construye para testear las hipótesis:
H0:β1=β2=···=βp−1=0
H1: no todos los βk (k=1,2,…,p−1) son iguales a 0.
Observemos que H0 dice que no hay vínculo entre la variable respuesta y las regresoras. En cambio, H1 dice que al menos una de las variables regresoras sirve para predecir a Y. Los resultados de este test se pueden observar haciendo un summary() del modelo o glance().
glance(modelo_ln_base)
F-statistic: 770.4 on 5 and 7018 DF, p-value: < 2.2e-16
Test significativo por lo que determinamos que se rechaza la hipotesis nula, nuestro modelo es relevante en comparación al modelo base, para determinar que arroja resultados significativos tomando la variable predictora
El R-cuadrado permite medir el porcentaje de variabilidad del fenómeno que el modelo logra explicar. Por este motivo es una métrica que nos permite evaluar la capacidad explicativa del modelo y poder comparar modelos entre sí bajo ciertas condiciones.
En nuestro caso tenemos un valor de Multiple R-squared: 0.3544, Adjusted R-squared: 0.3539
encuesta_salud_train %>% glimpse
Rows: 7,024
Columns: 17
$ record <int> 502, 26488, 31473, 14154, 36578, 53730, …
$ edad <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero <chr> "Femenino", "Masculino", "Masculino", "M…
$ nivel_educativo <chr> "2do año/11vo grado nivel Polimodal o 4t…
$ altura <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ frecuencia_hambre_mensual <chr> "Rara vez", "Rara vez", "Nunca", "Nunca"…
$ dias_consumo_comida_rapida <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1…
$ edad_consumo_alcohol <chr> "14 o 15 años", "7 años o menos", "Nunca…
$ consumo_diario_alcohol <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…
$ consumo_semanal_frutas <chr> "No comí frutas durante los últimos 7 dí…
$ consumo_semanal_verdura <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_gaseosas <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_snacks <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_comida_grasa <chr> "No comí comida alta en grasa en los últ…
$ id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1…
encuesta_salud_second_model = encuesta_salud_train %>%
# Seleccionamos las variables de interés
dplyr::select(edad, genero, altura, peso,consumo_semanal_snacks)
encuesta_salud_second_model %>% glimpse
Rows: 7,024
Columns: 5
$ edad <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, 16, 14,…
$ genero <chr> "Femenino", "Masculino", "Masculino", "Masculin…
$ altura <int> 165, 178, 172, 170, 170, 178, 156, 163, 164, 16…
$ peso <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, 100, 33…
$ consumo_semanal_snacks <chr> "1 a 3 veces durante los últimos 7 días", "No c…
Dejamos al valor No comí comida salada o snacks en los últimos 7 días de consumo semanal de snacks como categoría basal en una nueva variable que llamaremos consumo_semanal_snacks_relevel
encuesta_salud_second_model$consumo_semanal_snacks_relevel <- relevel(as.factor(encuesta_salud_second_model$consumo_semanal_snacks), ref = 8)
Graficámos la variable categórica consumoSemanalSnacks para analizar sus valores en función del peso
# armo boxplots paralelos consumo de snacks en función del peso
ggplot( encuesta_salud_second_model, aes(x = fct_reorder(consumo_semanal_snacks_relevel, peso, .desc = T), y = peso)) +
geom_boxplot(alpha = 0.75, aes(fill = consumo_semanal_snacks_relevel)) +
theme_minimal() +
theme(legend.position = 'none')+
labs(y = "Peso", x = "Consumo snacks") +
ggtitle("Boxplots del peso en función del consumo de snacks")+
theme (axis.text.x = element_text(face="italic", colour="dark grey", size = 8, angle = 90))
modelo_ln_2_r <- lm(peso ~ consumo_semanal_snacks_relevel + edad + genero + altura + (edad * genero), data = encuesta_salud_second_model)
# Resumen del modelo
tidy_ln_2_r <- tidy(modelo_ln_2_r, conf.int = TRUE)
print(tidy_ln_2_r)
# A tibble: 12 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -64.2 2.83 -22.7 4.08e-110 -69.7 -58.7
2 consumo_semanal_sn… -1.35 0.276 -4.89 1.03e- 6 -1.89 -0.810
3 consumo_semanal_sn… -0.608 0.456 -1.33 1.83e- 1 -1.50 0.286
4 consumo_semanal_sn… -1.09 0.685 -1.60 1.10e- 1 -2.44 0.247
5 consumo_semanal_sn… -1.28 1.01 -1.26 2.07e- 1 -3.26 0.707
6 consumo_semanal_sn… -2.27 0.450 -5.05 4.64e- 7 -3.15 -1.39
7 consumo_semanal_sn… -2.57 0.881 -2.91 3.59e- 3 -4.29 -0.840
8 consumo_semanal_sn… -4.44 1.95 -2.28 2.29e- 2 -8.27 -0.615
9 edad 1.22 0.121 10.1 7.55e- 24 0.986 1.46
10 generoMasculino -4.61 2.68 -1.72 8.52e- 2 -9.86 0.640
11 altura 0.643 0.0145 44.3 0 0.614 0.671
12 edad:generoMasculi… 0.391 0.179 2.19 2.88e- 2 0.0405 0.742
Intercept
Coeficiente edad:generoMasculino:
Coeficiente del consumo semanal de snacks:
El coeficiente estimado de consumo semanal snacks 1 a 3 veces durante los últimos 7 días es de -1.35, si mantenemos el resto de las variables constantes, el incremento de una unidad para consumo_semanal_snacks 1 a 3 veces durante los últimos 7 días, corresponde a una baja de 1.35kg, en promedio en el peso de la persona. => Significativa
El coeficiente estimado de consumo_semanal_snacks 1 vez al día es de -0.608, si mantenemos el resto de las variables constantes, el incremento de una unidad para consumo_semanal_snacks 1 vez al día corresponde a una baja de 0.608, en promedio en el peso de la persona. No significativa
El coeficiente estimado de consumo_semanal_snacks 2 veces al día es de -1.09, si mantenemos el resto de las variables constantes, el incremento de una unidad para consumo_semanal_snacks 2 veces al día corresponde a una baja de 1.09, en promedio en el peso de la persona. => No significativa
El coeficiente estimado de consumo_semanal_snacks 3 veces al día es de -1.28, si mantenemos el resto de las variables constantes, el incremento de una unidad para consumo_semanal_snacks 3 veces al día corresponde a una baja de 1.28, en promedio en el peso de la persona. => No significativa
El coeficiente estimado de consumo_semanal_snacks 4 a 6 veces durante los últimos 7 días es de -2.27, si mantenemos el resto de las variables constantes, el incremento de una unidad para consumo_semanal_snacks 4 a 6 veces durante los últimos 7 días corresponde a una baja de 2.27kg, en promedio en el peso de la persona. => Significativa
El coeficiente estimado de consumo_semanal_snacks 4 o más veces al día es de -2.57, si mantenemos el resto de las variables constantes, el incremento de una unidad para estimado de consumo_semanal_snacks 4 o más veces al día corresponde a una baja de 2.57kg, en promedio en el peso de la persona. => Significativa
El coeficiente estimado de consumo_semanal_snacks Dato perdido es de -4.44, si mantenemos el resto de las variables constantes, el incremento de una unidad para consumo_semanal_snacksDato perdido corresponde a una baja de 4.44kg, en promedio en el peso de la persona. => Significativa
Los coeficientes para las diferentes categorias de consumo semanal de snacks algunos son significativos y otros no.
Realizamos el test de anova para ver la significancia de las variables en su conjunto
tidy(anova(modelo_ln_2_r))
La tabla de ANOVA muestra que, según el resultado del test F, la variable consumo_semanal_snacks en su conjunto resulta estadísticamente significativa para explicar al precio (p-valor < 0.05).
Creamos una variable nueva que agrupe las tres variables categóricas que nos dieron una NO significancia en el análisis previo, las mismas serian
[3] “consumo_semanal_snacks_relevel1 vez al día”
[4] “consumo_semanal_snacks_relevel2 veces al día”
[5] “consumo_semanal_snacks_relevel3 veces al día”
encuesta_salud_train %>% glimpse
Rows: 7,024
Columns: 17
$ record <int> 502, 26488, 31473, 14154, 36578, 53730, …
$ edad <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero <chr> "Femenino", "Masculino", "Masculino", "M…
$ nivel_educativo <chr> "2do año/11vo grado nivel Polimodal o 4t…
$ altura <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ frecuencia_hambre_mensual <chr> "Rara vez", "Rara vez", "Nunca", "Nunca"…
$ dias_consumo_comida_rapida <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1…
$ edad_consumo_alcohol <chr> "14 o 15 años", "7 años o menos", "Nunca…
$ consumo_diario_alcohol <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…
$ consumo_semanal_frutas <chr> "No comí frutas durante los últimos 7 dí…
$ consumo_semanal_verdura <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_gaseosas <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_snacks <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_comida_grasa <chr> "No comí comida alta en grasa en los últ…
$ id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1…
encuesta_salud_second_model_agrupado = encuesta_salud_train %>%
# Seleccionamos las variables de interés
dplyr::select(edad, genero, altura, peso,consumo_semanal_snacks)
encuesta_salud_second_model_agrupado$consumo_semanal_snacks_relevel <- relevel(as.factor(encuesta_salud_second_model_agrupado$consumo_semanal_snacks), ref = 8)
Redefinimos la variable agrupando los valores categorizados como 1 vez al día, 2 veces al día y tres veces al día, en 1 a 3 veces al dia.
encuesta_salud_second_model_agrupado = encuesta_salud_second_model_agrupado %>% mutate(consumo_semanal_snacks = case_when(startsWith(consumo_semanal_snacks, "1 vez al día") ~ "1 a 3 veces al dia",
startsWith(consumo_semanal_snacks, "2 veces al día") ~ "consumo_semanal_snacks 1 a 3 veces al dia",
startsWith(consumo_semanal_snacks, "3 veces al día") ~ "consumo_semanal_snacks 1 a 3 veces al dia",
TRUE ~ consumo_semanal_snacks))
#encuesta_salud_second_model_agrupado$consumo_semanal_snacks
modelo_ln_2_r_agrupado <- lm(peso ~ consumo_semanal_snacks + edad + genero + altura + (edad * genero), data = encuesta_salud_second_model_agrupado)
# Resumen del modelo
tidy_ln_2_r_agrupado <- tidy(modelo_ln_2_r_agrupado, conf.int = TRUE)
print(tidy_ln_2_r_agrupado)
# A tibble: 11 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -64.8 2.83 -22.9 4.14e-112 -70.4 -59.3
2 consumo_semanal_sn… -0.744 0.440 -1.69 9.11e- 2 -1.61 0.119
3 consumo_semanal_sn… -1.66 0.565 -2.94 3.28e- 3 -2.77 -0.554
4 consumo_semanal_sn… -1.96 0.944 -2.07 3.81e- 2 -3.81 -0.108
5 consumo_semanal_sn… -0.542 0.676 -0.801 4.23e- 1 -1.87 0.784
6 consumo_semanal_sn… -3.83 1.98 -1.93 5.30e- 2 -7.72 0.0503
7 consumo_semanal_sn… 0.608 0.456 1.33 1.83e- 1 -0.286 1.50
8 edad 1.22 0.121 10.1 7.57e- 24 0.986 1.46
9 generoMasculino -4.61 2.68 -1.72 8.52e- 2 -9.86 0.640
10 altura 0.643 0.0145 44.3 0 0.614 0.671
11 edad:generoMasculi… 0.391 0.179 2.19 2.88e- 2 0.0405 0.742
Observamos que la redefinición de las variables no significativas generando un agrupamiento para de las mismas mantiene la no significancia de estos valores para predecir el peso de una persona.
Ejecutamos las funcines tidy, glance y summary para obtener tanto la significacia de los atributos del modelo, la significancia del modelo en si (Test F) y la variabilidad del mismo (R-cuadrado y R-cuadrado ajustado).
options("scipen"=1)
tidy_ln_2_r %>%
dplyr::select(term, statistic, p.value, conf.low, conf.high)
glance(modelo_ln_2_r)
summary(modelo_ln_2_r)
Call:
lm(formula = peso ~ consumo_semanal_snacks_relevel + edad + genero +
altura + (edad * genero), data = encuesta_salud_second_model)
Residuals:
Min 1Q Median 3Q Max
-29.005 -6.452 -1.157 5.052 75.767
Coefficients:
Estimate
(Intercept) -64.19940
consumo_semanal_snacks_relevel1 a 3 veces durante los últimos 7 días -1.35163
consumo_semanal_snacks_relevel1 vez al día -0.60788
consumo_semanal_snacks_relevel2 veces al día -1.09466
consumo_semanal_snacks_relevel3 veces al día -1.27596
consumo_semanal_snacks_relevel4 a 6 veces durante los últimos 7 días -2.27004
consumo_semanal_snacks_relevel4 o más veces al día -2.56697
consumo_semanal_snacks_relevelDato perdido -4.44042
edad 1.22309
generoMasculino -4.61150
altura 0.64292
edad:generoMasculino 0.39146
Std. Error
(Intercept) 2.82849
consumo_semanal_snacks_relevel1 a 3 veces durante los últimos 7 días 0.27645
consumo_semanal_snacks_relevel1 vez al día 0.45601
consumo_semanal_snacks_relevel2 veces al día 0.68457
consumo_semanal_snacks_relevel3 veces al día 1.01150
consumo_semanal_snacks_relevel4 a 6 veces durante los últimos 7 días 0.44992
consumo_semanal_snacks_relevel4 o más veces al día 0.88112
consumo_semanal_snacks_relevelDato perdido 1.95125
edad 0.12102
generoMasculino 2.67877
altura 0.01453
edad:generoMasculino 0.17902
t value
(Intercept) -22.697
consumo_semanal_snacks_relevel1 a 3 veces durante los últimos 7 días -4.889
consumo_semanal_snacks_relevel1 vez al día -1.333
consumo_semanal_snacks_relevel2 veces al día -1.599
consumo_semanal_snacks_relevel3 veces al día -1.261
consumo_semanal_snacks_relevel4 a 6 veces durante los últimos 7 días -5.045
consumo_semanal_snacks_relevel4 o más veces al día -2.913
consumo_semanal_snacks_relevelDato perdido -2.276
edad 10.106
generoMasculino -1.722
altura 44.250
edad:generoMasculino 2.187
Pr(>|t|)
(Intercept) < 2e-16
consumo_semanal_snacks_relevel1 a 3 veces durante los últimos 7 días 1.03e-06
consumo_semanal_snacks_relevel1 vez al día 0.18257
consumo_semanal_snacks_relevel2 veces al día 0.10985
consumo_semanal_snacks_relevel3 veces al día 0.20719
consumo_semanal_snacks_relevel4 a 6 veces durante los últimos 7 días 4.64e-07
consumo_semanal_snacks_relevel4 o más veces al día 0.00359
consumo_semanal_snacks_relevelDato perdido 0.02289
edad < 2e-16
generoMasculino 0.08520
altura < 2e-16
edad:generoMasculino 0.02880
(Intercept) ***
consumo_semanal_snacks_relevel1 a 3 veces durante los últimos 7 días ***
consumo_semanal_snacks_relevel1 vez al día
consumo_semanal_snacks_relevel2 veces al día
consumo_semanal_snacks_relevel3 veces al día
consumo_semanal_snacks_relevel4 a 6 veces durante los últimos 7 días ***
consumo_semanal_snacks_relevel4 o más veces al día **
consumo_semanal_snacks_relevelDato perdido *
edad ***
generoMasculino .
altura ***
edad:generoMasculino *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.887 on 7012 degrees of freedom
Multiple R-squared: 0.3585, Adjusted R-squared: 0.3575
F-statistic: 356.3 on 11 and 7012 DF, p-value: < 2.2e-16
El porcentaje de variabilidad que explica este nuevo modelo es de R-cuadrado = 0.359, R-cuadrado ajustado = 0.358 Evaluando la significativida del modelo, obtenemos F-statistic: 356.3 on 11 and 7012 DF, p-value: < 2.2e-16 Modelo significativo
Proponemos un modelo que evalue el peso a partir del nivel educativo, consumo_semanal_comida_grasa, edad, genero, altura y una interrelación de altura, edad y género.
f(peso) ~ nivel_educativo + consumo_semanal_comida_grasa + edad + genero + altura + (altura * edad * genero)
Seleccionamos nuestras variables de interés
encuesta_salud_train %>% glimpse
Rows: 7,024
Columns: 17
$ record <int> 502, 26488, 31473, 14154, 36578, 53730, …
$ edad <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero <chr> "Femenino", "Masculino", "Masculino", "M…
$ nivel_educativo <chr> "2do año/11vo grado nivel Polimodal o 4t…
$ altura <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ frecuencia_hambre_mensual <chr> "Rara vez", "Rara vez", "Nunca", "Nunca"…
$ dias_consumo_comida_rapida <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1…
$ edad_consumo_alcohol <chr> "14 o 15 años", "7 años o menos", "Nunca…
$ consumo_diario_alcohol <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…
$ consumo_semanal_frutas <chr> "No comí frutas durante los últimos 7 dí…
$ consumo_semanal_verdura <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_gaseosas <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_snacks <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_comida_grasa <chr> "No comí comida alta en grasa en los últ…
$ id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1…
encuesta_salud_third_model = encuesta_salud_train %>%
# Seleccionamos las variables de interés
dplyr::select(edad, genero, altura, peso,nivel_educativo,consumo_semanal_comida_grasa)
encuesta_salud_third_model %>% glimpse
Rows: 7,024
Columns: 6
$ edad <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, 1…
$ genero <chr> "Femenino", "Masculino", "Masculino", "Ma…
$ altura <int> 165, 178, 172, 170, 170, 178, 156, 163, 1…
$ peso <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, 1…
$ nivel_educativo <chr> "2do año/11vo grado nivel Polimodal o 4to…
$ consumo_semanal_comida_grasa <chr> "No comí comida alta en grasa en los últi…
Generamos el modelo lineal
modelo_ln_3_r <- lm(peso ~ nivel_educativo + consumo_semanal_comida_grasa + edad + genero + altura + (altura * edad * genero), data = encuesta_salud_third_model)
# Resumen del modelo
tidy_ln_3_r <- tidy(modelo_ln_3_r, conf.int = TRUE)
print(tidy_ln_3_r)
# A tibble: 20 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -32.9 38.6 -0.852 0.394 -1.09e+2 42.8
2 nivel_educativo2do a… -0.779 0.356 -2.19 0.0287 -1.48e+0 -0.0808
3 nivel_educativo3er a… -0.423 0.411 -1.03 0.304 -1.23e+0 0.383
4 nivel_educativo8vo g… -0.930 0.458 -2.03 0.0424 -1.83e+0 -0.0321
5 nivel_educativo9no g… -0.400 0.379 -1.06 0.290 -1.14e+0 0.342
6 nivel_educativoDato … -1.38 1.02 -1.35 0.176 -3.39e+0 0.622
7 consumo_semanal_comi… 1.10 0.428 2.57 0.0102 2.61e-1 1.94
8 consumo_semanal_comi… -1.42 0.630 -2.26 0.0239 -2.66e+0 -0.188
9 consumo_semanal_comi… -1.33 0.986 -1.35 0.177 -3.27e+0 0.601
10 consumo_semanal_comi… -0.915 0.349 -2.62 0.00870 -1.60e+0 -0.231
11 consumo_semanal_comi… -0.0676 0.849 -0.0796 0.937 -1.73e+0 1.60
12 consumo_semanal_comi… -2.94 1.40 -2.10 0.0354 -5.68e+0 -0.200
13 consumo_semanal_comi… 0.0307 0.308 0.0996 0.921 -5.74e-1 0.635
14 edad -0.343 2.57 -0.133 0.894 -5.38e+0 4.70
15 generoMasculino 11.1 50.3 0.221 0.825 -8.75e+1 110.
16 altura 0.440 0.241 1.82 0.0682 -3.29e-2 0.913
17 edad:altura 0.0101 0.0161 0.629 0.529 -2.14e-2 0.0416
18 generoMasculino:altu… -0.0673 0.308 -0.219 0.827 -6.71e-1 0.536
19 edad:generoMasculino -1.74 3.38 -0.516 0.606 -8.37e+0 4.88
20 edad:generoMasculino… 0.0111 0.0206 0.538 0.591 -2.93e-2 0.0515
glance(modelo_ln_3_r)
tidy(anova(modelo_ln_3_r))
Proponemos un modelo que evalue el peso a partir del nivel educativo, consumo_semanal_comida_grasa, genero, altura y una interrelación de dias consumo comida rápida y genero f(peso) ~ nivel_educativo + consumo_semanal_comida_grasa + genero + altura + (dias_consumo_comida_rapida * genero)
Seleccionamos nuestras variables de interés
encuesta_salud_train %>% glimpse
Rows: 7,024
Columns: 17
$ record <int> 502, 26488, 31473, 14154, 36578, 53730, …
$ edad <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero <chr> "Femenino", "Masculino", "Masculino", "M…
$ nivel_educativo <chr> "2do año/11vo grado nivel Polimodal o 4t…
$ altura <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ frecuencia_hambre_mensual <chr> "Rara vez", "Rara vez", "Nunca", "Nunca"…
$ dias_consumo_comida_rapida <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1…
$ edad_consumo_alcohol <chr> "14 o 15 años", "7 años o menos", "Nunca…
$ consumo_diario_alcohol <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…
$ consumo_semanal_frutas <chr> "No comí frutas durante los últimos 7 dí…
$ consumo_semanal_verdura <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_gaseosas <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_snacks <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_comida_grasa <chr> "No comí comida alta en grasa en los últ…
$ id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1…
encuesta_salud_fourth_model = encuesta_salud_train %>%
# Seleccionamos las variables de interés
dplyr::select(edad, genero, altura, peso,nivel_educativo,dias_consumo_comida_rapida)
encuesta_salud_fourth_model %>% glimpse
Rows: 7,024
Columns: 6
$ edad <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, 16,…
$ genero <chr> "Femenino", "Masculino", "Masculino", "Masc…
$ altura <int> 165, 178, 172, 170, 170, 178, 156, 163, 164…
$ peso <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, 100…
$ nivel_educativo <chr> "2do año/11vo grado nivel Polimodal o 4to a…
$ dias_consumo_comida_rapida <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1, 3…
Generamos nuestro modelo lineal
modelo_ln_4_r <- lm(peso ~ nivel_educativo + dias_consumo_comida_rapida + genero + altura + (dias_consumo_comida_rapida * genero), data = encuesta_salud_fourth_model)
# Resumen del modelo
tidy_ln_4_r <- tidy(modelo_ln_4_r, conf.int = TRUE)
print(tidy_ln_4_r)
# A tibble: 10 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -5.12e+1 2.31 -22.2 3.15e-105 -55.7 -46.6
2 nivel_educativo2do… 1.86e-1 0.348 0.534 5.93e- 1 -0.496 0.867
3 nivel_educativo3er… 1.50e+0 0.372 4.03 5.52e- 5 0.771 2.23
4 nivel_educativo8vo… -2.76e+0 0.425 -6.48 1.00e- 10 -3.59 -1.92
5 nivel_educativo9no… -1.52e+0 0.370 -4.10 4.20e- 5 -2.24 -0.791
6 nivel_educativoDat… -8.26e-1 1.03 -0.801 4.23e- 1 -2.85 1.20
7 dias_consumo_comid… -5.82e-3 0.117 -0.0498 9.60e- 1 -0.235 0.223
8 generoMasculino 1.41e+0 0.306 4.62 3.86e- 6 0.815 2.01
9 altura 6.72e-1 0.0143 47.0 0 0.644 0.700
10 dias_consumo_comid… -4.14e-1 0.175 -2.37 1.79e- 2 -0.757 -0.0712
glance(modelo_ln_4_r)
tidy(anova(modelo_ln_4_r))
Agrupamos todos los modelos que nos interesa comparar, en este caso es el modelo inicial, el modelo categóricas con las categorías redefinidas de la variable consumoSemanalSnacks y los dos modelos desarrollados en este punto, interrelación (altura * edad * genero) y (dias_consumo_comida_rapida * genero)
# armamos lista con todos los modelos
models <- list(modelo_ln_4_r = modelo_ln_4_r, modelo_ln_3_r = modelo_ln_3_r, modelo_ln_2_r = modelo_ln_2_r, modelo_ln_base = modelo_ln_base)
# calculamos las variables resumen
map_df(models, tidy, .id = "model")
Utilizamos la función augment para predecir el peso sobre el dataset de training
lista_predicciones_training = map(.x = models, .f = augment)
Calculamos el RMSE para todos los modelos a evaluar
map_dfr(.x = lista_predicciones_training, .f = rmse, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)
Calculamos el MAE para todos los modelos a evaluar
map_dfr(.x = lista_predicciones_training, .f = mae, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)
Obtenemos el R-cuadrado y R-cuadrado ajustado para poder comparar la variabilidad de los modelos
df_evaluacion_train = map_df(models, glance, .id = "model") %>%
# ordenamos por R2 ajustado
arrange(desc(adj.r.squared))
df_evaluacion_train
Concluimos sobre el dataset de training
Utilizamos la función augment para predecir el peso sobre el dataset de training
# Aplicamos la función augment a los 4 modelos con el set de testing
encuesta_salud_test$consumo_semanal_snacks_relevel <- relevel(as.factor(encuesta_salud_test$consumo_semanal_snacks), ref = 8)
lista_predicciones_testing = map(.x = models, .f = augment, newdata = encuesta_salud_test)
Calculamos el RMSE para todos los modelos a evaluar con las predicciones en training
# Obtenemos el RMSE para los 4 modelos
map_dfr(.x = lista_predicciones_testing, .f = rmse, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)
Calculamos el MAE para todos los modelos a evaluar con las predicciones en training
map_dfr(.x = lista_predicciones_testing, .f = mae, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)
Concluimos sobre el dataset de testing:
Realizaremos diferentes pruebas para validar el cumplimiento o no de los supuestos del modelo lineal, estos consisten de: εi∼N(0,σ2) independientes entre sí. Los errores son inobservables, por lo tanto tendremos que trabajar con su correlato empírico: los residuos en las técnicas de diagnóstico.
Procedemos a graficar el comportamiento de los residuos del modelo
au_modelos = map_df(models, augment, .id = "model")
g1 = ggplot(au_modelos %>% filter(model == "modelo_ln_base"),
aes(.fitted, .resid)) +
geom_point() +
geom_hline(yintercept = 0) +
geom_smooth(se = FALSE) +
labs(title = "Residuos vs valores predichos MODELO BASE") +
theme_bw()
g2 = ggplot(au_modelos %>% filter(model == "modelo_ln_base"),
aes(sample = .std.resid)) +
stat_qq() +
geom_abline() +
labs(title = "Normal QQ plot MODELO BASE") +
theme_bw()
g3 = ggplot(au_modelos %>% filter(model == "modelo_ln_base"),
aes(.fitted, sqrt(abs(.std.resid)))) +
geom_point() +
geom_smooth(se = FALSE) +
theme_bw() +
labs(title = "Scale-location plot MODELO BASE")
g4 = ggplot(au_modelos %>% filter(model == "modelo_ln_base"),
aes(.hat, .std.resid)) +
geom_vline(size = 2, colour = "white", xintercept = 0) +
geom_hline(size = 2, colour = "white", yintercept = 0) +
geom_point() +
geom_smooth(se = FALSE) +
theme_bw() +
labs(title = "Residual vs leverage MODELO BASE")
# grafico todos juntos
grid.arrange(g1, g2, g3, g4, nrow = 2)
Conclusiones:
Diagnóstico del modelo: El modelo creado no cumple con los supuestos del modelo lineal. Parecen existir problemas de falta de normalidad y presencia de observaciones de alto leverage.
Leemos el nuevo archivo con la incorporación de algunas observaciones adicionales que pueden incluir valores atípicos
library(MASS)
encuesta_salud_train_original = read_csv("/Users/dfontenla/Maestria/2022C2/EEA/practica/repo/EEA-2022/TP1/encuesta_salud_modelo6.csv") %>% mutate(id = 1:nrow(.))
encuesta_salud_train_original %>% glimpse()
Rows: 7,129
Columns: 17
$ record <dbl> 48568, 39818, 11360, 50320, 26325, 44811…
$ edad <dbl> 17, 14, 16, 14, 15, 13, 16, 13, 15, 14, …
$ genero <chr> "Femenino", "Masculino", "Masculino", "M…
$ nivel_educativo <chr> "3er año/12vo grado nivel Polimodal o 5t…
$ altura <dbl> 160, 189, 166, 178, 167, 164, 162, 170, …
$ peso <dbl> 51, 75, 50, 80, 50, 55, 60, 62, 48, 71, …
$ frecuencia_hambre_mensual <chr> "Nunca", "Nunca", "Nunca", "Nunca", "Nun…
$ dias_consumo_comida_rapida <dbl> 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0…
$ edad_consumo_alcohol <chr> "12 o 13 años", "Nunca tomé alcohol más …
$ consumo_diario_alcohol <dbl> 5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 0.0, …
$ dias_actividad_fisica_semanal <dbl> 2, 0, 0, 2, 1, 3, 0, 2, 3, 7, 2, 2, 1, 6…
$ consumo_semanal_frutas <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_verdura <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_gaseosas <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_snacks <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_comida_grasa <chr> "4 a 6 veces durante los últimos 7 días"…
$ id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1…
Graficamos la relación peso / altura para observar las anomalias insertadas en el nuevo dataset
encuesta_salud_train_original %>% ggplot(., aes(x = altura, y = peso)) +
geom_point(alpha=0.5) + #capa de los datos
theme_bw() +
labs(title="Modelo Lineal Múltiple: Altura x Peso", x="Altura", y="Peso")
Observamos que existen una cantidad de puntos relevantes que se alejan de forma categórica de la distribución general, una importante cantidad de outliers
Nos quedamos con las variables de interés para reproducir el modelo base
encuesta_salud_train_original_base = encuesta_salud_train_original %>%
# Seleccionamos las variables de interés
dplyr::select(edad, genero, altura, peso,dias_actividad_fisica_semanal,consumo_diario_alcohol)
encuesta_salud_train_original_base %>% glimpse
Rows: 7,129
Columns: 6
$ edad <dbl> 17, 14, 16, 14, 15, 13, 16, 13, 15, 14, …
$ genero <chr> "Femenino", "Masculino", "Masculino", "M…
$ altura <dbl> 160, 189, 166, 178, 167, 164, 162, 170, …
$ peso <dbl> 51, 75, 50, 80, 50, 55, 60, 62, 48, 71, …
$ dias_actividad_fisica_semanal <dbl> 2, 0, 0, 2, 1, 3, 0, 2, 3, 7, 2, 2, 1, 6…
$ consumo_diario_alcohol <dbl> 5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 0.0, …
Graficamos las relaciones entre variables del nuevo dataset
encuesta_salud_train_original_base %>% ggpairs(aes(color=genero), upper = list(continuous = wrap("cor", size = 3, hjust=0.8, alignPercent=0.15)), legend = 25) +
theme_bw() +
theme(axis.text.x = element_text(angle=45, vjust=0.5), legend.position = "bottom")
Representamos la correlación de Pearson entre variables
encuesta_salud_train_original_base %>%
correlate() %>% # convierte la matriz de corr en dataframe
shave() %>% # solo muestra información debajo de la diagonal principal
fashion() # acomoda los datos en forma tidy (por ej. redondeo de decimales)
Podemos observar que la correlación entre peso y altura en relación al dataset base, disminuyó considerablemente pasando de un 0.5 a un 0.29, gráficamente se visualiza una clara correspondencia entre estos valores, ahora bien para el nuevo dataset existe una cantidad de datos outliers significativos
Ajustamos el dataset a nuestro modelo lineal múltiple base (peso ~ edad + genero + altura + dias_actividad_fisica_semanal + consumo_diario_alcohol)
# ajustamos modelo lineal múltiple
modelo_ds_original_ln_base <- lm(peso ~ edad + genero + altura + dias_actividad_fisica_semanal + consumo_diario_alcohol, data = encuesta_salud_train_original_base)
# Resumen del modelo
tidy_ln_ds_original_base <- tidy(modelo_ds_original_ln_base, conf.int = TRUE)
tidy_ln_ds_original_base
Realizamos un diagnóstico del modelo a través de los residuos
plot(modelo_ds_original_ln_base)
Vemos la significancia de los atributos
options("scipen"=1)
tidy_ln_ds_original_base %>%
dplyr::select(term, statistic, p.value, conf.low, conf.high)
Analizamos la variabilidad explicativa del modelo
glance(modelo_ds_original_ln_base)
Observamos que la variabilidad porcentual explicada por R-cuadrado y R-cuadrado ajustado disminuye considerablemente en relación al mismo modelo entrenado en el dataset sin outliers
tidy(anova(modelo_ds_original_ln_base))
Agrupamos los modelos para poder compararlos
# armamos lista con todos los modelos
models_leverage <- list(modelo_ds_original_ln_base = modelo_ds_original_ln_base, modelo_ln_base = modelo_ln_base)
# calculamos las variables resumen
map_df(models_leverage, tidy, .id = "model")
lista_predicciones_training_leverage = map(.x = models_leverage, .f = augment)
Utilizamos la función glance sobre todos los modelos para comparar la variabilidad de los mismos
df_evaluacion_train_leverage = map_df(models_leverage, glance, .id = "model") %>%
# ordenamos por R2 ajustado
arrange(desc(adj.r.squared))
df_evaluacion_train_leverage
Calculamos el RMSE para todos los modelos a evaluar
map_dfr(.x = lista_predicciones_training_leverage, .f = rmse, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)
Calculamos el MAE para todos los modelos a evaluar
map_dfr(.x = lista_predicciones_training_leverage, .f = mae, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)
Realizando la comparación de las metricas R-cuadrado ajustado, RMSE y MAE, observamos que:
Concluimos con lo observado que nuestro modelo es sensible a outliers, y ante la precencia de estos, nuestras métricas de performance del modelo empeoran drásticamente
Ante la apreciación de como empeora nuestro modelo ante la existencia de outliers en nuestro dataset, procedemos a realizar un entrenamiento con un modelo lineal robusto el cual esperaremos que brinde mejores resultados debido a la ponderación que realizan sobre la influencia de los casos atípicos
# ajustamos modelo lineal múltiple
modelo_ds_original_ln_robust_base <- rlm(peso ~ edad + genero + altura + dias_actividad_fisica_semanal + consumo_diario_alcohol, data = encuesta_salud_train_original_base)
# # Resumen del modelo
tidy(modelo_ds_original_ln_base)
tidy(modelo_ds_original_ln_robust_base)
Comparando coeficientes, podemos ver como los relacionados al modelo robusto, no sufren grandes modificaciones en relación al modelo entrenado con el dataset sin outliers, mientras que en el modelo lineal común entrenado con el dataset con outliers, los coeficientes cambian considerablemente debido a estos datos.
options("scipen"=1)
tidy_ln_ds_original_robust <- tidy(modelo_ds_original_ln_robust_base, conf.int = TRUE)
tidy_ln_ds_original_robust %>% dplyr::select(term, statistic, conf.low, conf.high)
models_robust <- list(modelo_ln_5_robust = modelo_ds_original_ln_robust_base, modelo_ln_5_base = modelo_ds_original_ln_base)
robust_predictions_over_training = map(.x = models_robust, .f = augment)
map_dfr(.x = robust_predictions_over_training, .f = rmse, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)
map_dfr(.x = robust_predictions_over_training, .f = mae, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)
Evaluamos la performance de los modelos sobre el conjunto de testing
# Aplicamos la función augment a los 4 modelos con el set de testing
lista_predicciones_testing_robust = map(.x = models_robust, .f = augment, newdata = encuesta_salud_test)
# Obtenemos el RMSE para los 4 modelos
map_dfr(.x = lista_predicciones_testing_robust, .f = rmse, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)
map_dfr(.x = lista_predicciones_testing_robust, .f = mae, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)
Evaluando los modelos sobre el conjunto de pruebas
Conclusión, ante la presencia de una cantidad significativa de outliers, utilizar un modelo lineal robusto nos sirvio para poder tener un modelo que ante la presencia de estos valores atípicos se sigue comportando de una manera similar al modelo que fue entrenado sin outliers, mientras que el modelo lineal no robusto, siendo entrenado en un dataset con outliers, cambia significativamente sus coeficientes. Hace sentido poner foco sobre la métria MAE ya que es más robusto cuando los datos tienen outliers o datos atípicos y es la mejor opción a usar en esos casos